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ABSTRACT 

The nature of vertical heating of disk stars in the inner as well as the outer region of disk galaxies is 
studied. The galactic bar (which is the strongest non-axisymmetric pattern in the disk) is shown to be 
a potential source of vertical heating of the disk stars in the inner region. Using a nearly self-consistent 
high-resolution N-body simulation of disk galaxies, the growth rate of the bar potential is found to 
be positively correlated with the vertical heating exponent in the inner region of galaxies. We also 
characterize the vertical heating in the outer region where the disk dynamics is often dominated by 
the presence of transient spiral waves and mild bending waves. Our simulation results suggest that 
the non-axisymmetric structures are capable of producing the anisotropic heating of the disk stars. 

Subject headings: galaxies: evolution - galaxies: kinematics and dynamics - galaxies: structure - 
galaxies: spiral - stellar dynamics - galaxies: haloes 



1. INTRODUCTION 

The random velocities of disk stars in the solar neigh- 
borhood are known to increase with time since the 
1940's and this fact has been verified by recent Hippar- 

cos observation s. The analysis of the Hipparcos data 

(jBinnev et al. 1 120001 ) reveals that the stellar velocity 
dispersions increas e continuously (rather th an episodi- 
cally as claimed bv lEdvardsson et al. I (jl993f )) with time 
which is the well-known disk heating problem in our 
Galaxy. The phenomenon of disk heating is not spe- 
cific for our Galaxy since the heating of stellar disks 
is also known to persist in external gal axies such as in 
NGC 2985, NGC 2460, an d NGC 2775 (IGerssen et al. I 
120001 iShapiro et al. 1 120031 ) . Although mechanisms have 
been identified for the disk heating, a full understanding 
has yet to be achieved. It is known that disk heating is 
anisotropic; that is, disk stars are preferentially heated 
along the plane in comparison to its perpendicular di- 
rection. Thus, the main problem in disk heating is two- 
fold. The first is the requirement to explain the contin- 
ual heating of disk stars especially in the case of isolated 
galaxies for which tidal heating or external perturbations 
can be ruled out. Secondly, such a mechanism must be 
anisotropic in nature. Taken together, disk heating re- 
mains a funda mental problem in galactic dynamics . 

Historically, iSpitzer fc Schwarschild 1 (|195lL 119531 ) first 
showed that the secular increase in the stellar velocity 
dispersions could arise as a result of the scattering of 
disk stars with giant molecular clouds (GMCs) based on 
two dimensional calculations. Subsequently, there has 
been considerable development and progress in under- 
standi ng the d i sk he ating problem in galaxies. For ex- 
ample, lLacev I (|1984D extended the previous theoretical 
work to a three dimensional disk and found that the 
velocity dispersions increase with time, t, as V where 
7 = 0.25. In addition, it w as demonstrated that tran- 
sient stochastic spiral arms (jBarbanis fc Woltier Ill967t 
ICarlberg fc Sellwood 1 119851 ) could also produce secular 
heating of the disk stars, however, such effects were only 



effective in the plane of the galaxy. The effect of gi- 
ant molecular clouds in combination with the effect of 
stochastic spirals could s catter the stars off the plane 
([Jenkins fc Binneyl 119901 ) to account for both the ra- 
dial and vertical heating. Alternatively, using an exten- 
sive orbit anal ysis of the s tars subject to a 3D galactic 
bar potential, [Pfcnnigcr (1984, 119851) has shown that 
a large fraction of the phase space is chaotic and sug- 
gested that all the disk stars trapped in the chaotic phase 
space would be heated up in all directions. Other possi- 
ble candidate heating mechanisms include massive dark 
halo objec ts, e.g. black holes or dark clusters of mass 
- 10 6 M (lLacev fc Ostlik^[l98irCarr fc Lacevll987l) 
or a combi nation of giant molecular clo uds and halo 
black holes ([Hanninen fc Flvn 112001 12004D . Within the 
framework of standard cosmological models, a spectrum 
of su bhalos and their substructures could heat up the 
disk (iSanchez-Sakedo I [T999t lArdi et al. I [20031 ) . How- 
ever, iFont et al. I ([20011 ) found that the subhalos are not 
efficient perturbers for disk heating because the orbits 
of the subhalos rarely take them near the disk. Finally, 
satellite infall onto the gal actic disk could pr oduce an 
abrupt heating of the disk ([Quinn et al. 1119931 ) or mod- 
erate heating could arise due to sinking satellites on pro- 
grade orbits ([Velazquez fc White IH999I ). The many pos- 
sible heat ing mechanisms o f the disk stars has been re- 
viewed bv lPfenniger I (|1993l ) . Note that while most mech- 
anisms produce radial heating of the disk stars they have 
not been shown to provide a satisfactory explanation of 
the detailed nature of vertical heating. In fact, studies fo- 
cusing on the vertical heating are comparatively lacking. 
To be more precise, a generic and robust internal source 
remains to be identified which could be responsible for 
the continual vertical heating of the disk. 

Here, we primarily concentrate on the issue of verti- 
cal disk heating. Galactic disks often consist of various 
non-axisymmetric patterns such as bars, spirals, warps, 
corrugations, and rings. These patterns may be grow- 
ing or transient and, in general, may be time-dependent. 
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It is natural to ask whether such patterns contribute to 
the disk vertical heating. That is, what is the nature 
of vertical heating due, for example, to a growing bar? 
How does a growing bar heat the disk stars in the verti- 
cal direction? Heating due to non-axisymmetric patterns 
are important for isolated galaxies as even for the Milky 
Way, there are prominent an d characteristic signatures of 
non-axisymmetric patterns (Dehncn 1998) in the phase 
space in the solar neighbourhood. 

In this paper, we explore possible mechanisms giving 
rise to continuous heating of the disk stars in the ver- 
tical direction. In particular, we show that stellar bars 
that are formed nearly spontaneously in disk galaxies 
are capable of efficiently scattering stars in the vertical 
direction. Since transient spiral arms are often associ- 
ated with bars, they also contribute to the disk heat- 
ing. In order to understand the nature of vertical heat- 
ing of the disk stars, we have performed a large number 
(statistically meaningful) of N-body simulations of model 
disk galaxies constructed over a wide range of parameter 
space. 

The paper is organized as follows. Section 2 describes 
the galaxy models used in the simulation. A description 
of the N-body simulation and the physical basis for the 
choice of our parameter space is presented in section 3. 
The heating model used to interpret the numerical simu- 
lations is described in section 4. Section 5 describes the 
heating due to non-axisymmetric structures. The results 
on the correlation studies are presented in section 6. The 
discussion and conclusions are presented in section 7 and 
8 respectively. 

2. GALAXY MODELS 

The construction of a very thin equilibrium disk model 
involves a non-trivial procedure in N-body simulation. 
H owever, using the self- c onsist ent bulge-disk-halo model 
of IKuiiken fc Dubinski I (fl995l hence KD95) we are able 
to construct extremely thin equilibrium model disks 
(with the ratio of scale-height to scale length ~ 0.01). 
Their prescription provides nearly exact solution of the 
collisionless Boltzmann and Poisson equations which are 
suitable for studying disk stability related problems, al- 
lowing one to construct a wide range of initial models 
from a large parameter space. All the components in our 
models are active (i.e., the gravitational potential of each 
component can respond to an external or internal pertur- 
bation) and, hence, provide a realistic representation for 
the evolution and structure of the galaxies. Below we 
briefly describe each component of the model separately 
for the sake of completeness. For more details, the reader 
is referred to KD95. 

A sp herical live bulge is constr ucted from the King 
model (jBinnev fc Tremaine 1 1 1987h and the correspond- 
ing distribution function (DF) is given by 

r p b (2™ fc 2 )- 3 /v*°-*^ 

fb(E)= ! x{e -(E-^ c )/al _ 1} HE<y c , (1) 

[ otherwise. 

Here, the bulge is simulated using three parameters, 
namely the cut-off potential (^ c ), central bulge density 
(pb) and en, governing the central bulge velocity disper- 
sion. The depth of the potential well is measured by ^q. 
An axisymmetric live dark matter halo is simulated 



using the distribution function of a lowered Evans model 
and is given as 

f [(AL 2 z + B)e- E /°l+C] 
f dm (E, L 2 Z ) = < x {e- E /°i - 1) he < 0, 

[_ otherwise. 

(2) 

The velocity and density scales are given by au and p\ 
respectively. The halo core radius R c and the flattening 
parameter q together with p\ are contained in the param- 
eters A,B, and C. All the simulated halos are oblate in 
shape and kept at a constant value q — 0.8 for simplicity. 

The disk distribution function is constructed using the 
approximate third integral given by E z — ^v^+^(R, z) — 
^(R, 0), the energy of the vertical oscillations. This third 
integral is approximately conserved for orbits near the 
disk mid-plane. The radial density of the disk is approx- 
imately exponential with a truncation and the vertical 
density is chosen to depend exponentially on the verti- 
cal potential * 2 (-R, z) = ^{R, z) - ty(R, 0). The volume 
density of the axisymmetric disk is given by 

p J R , z ) = _^* e -R/R* erfc ( = - ~ Rout ) f d (z), 

(3) 

where f d (z) = exp(-0.8676* 2 (i?, z)/V z (R, h z )) governs 
the vertical structure of the disk, erfc is the complemen- 
tary error function. In the above equation, M d is the disk 
mass, Rd is the scale length and h z is the scale height. 

There are 13 parameters required to construct a par- 
ticular galaxy model. The galactic disk produced in this 
manner remains in equilibrium as long as the simula- 
tion runs. The dark matter halo generated using the 
lowered Evans model (jEvans lfl993h has a constant den- 
sity core which is probably appropriate for the low sur- 
face brightness galaxies which we simulate. Each model 
galaxy is constructed such that an initially almost flat 
rotation curve is produced. The disk outer radius (R ou t) 
is fixed at about 6.5Rd and a truncation width ~ 0.3Rd 
is adopted within which the disk density smoothly de- 
creases to zero at the outer radius. 

3. N-BODY SIMULATION 

The primary aim of the present study is to achieve an 
understanding of the nature and origin of vertical heating 
of disk stars. In order to gain the necessary insight, we 
perform a large number of simulations of isolated galax- 
ies. The models are evolved for a sufficiently long period 
of time so as to examine the nature of vertical heating in 
different spatial regions of the disk. 

Using the KD95 method, we build the initial condi- 
tions for all the model galaxies. The method is not 
fully self-consistent because the disk distribution func- 
tion is not known exactly. Each model galaxy con- 
sists of a bulge, disk and dark matter halo as described 
above. Since the initial conditions are constructed based 
on distribution functions, they are suitable for study- 
ing the long-term evolution of the non-axisymmetric 
patterns in the disk. Building initial galaxy models 
with prescribed galaxy properties such as the Toomre 
Q(r) — oy (r)K(r)/3.36G£(r), ratio of dark-to-disk mass 
(Mh/Md) or ratio of velocity dispersions (cr z /cr r ) is not 
straightforward because the bulge and the halo models 
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are not derived using their mass profiles but from their 
distribution functions. Hence, the relation between the 
Toomre Q parameter and Mh /Md or between o~ z / ay and 
Mh/Md is not known from the outset. Here, a r ,a z ,n, 
and £ denote the velocity dispersion in the radial and 
vertical direction, epicyclic frequency, and surface den- 
sity of stars at a given radius respectively. 

The disk, bulge and halo parameters are chosen such 
that an almost flat rotation curve is produced in the 
outer region. We assume a scale length for the disk to be 
Rd = 4.0 kpc. The radial and vertical forces are normal- 
ized such that it produces the circular velocity V c = 200 
kms^ 1 at about 2 disk scale lengths. The unit of time 
is 2.1 x 10 7 yr. The initial disk scale height (h z ) varies 
from 40 pc to 400 pc in our sample of model galaxies. 
The softening length for the disk particles is 10 pc, for 
the bulge particles 40 pc and for the halo particles the 
softening length is 36 pc. A total of 2.2 x 10 6 to 1.0 x 10 7 
particles has been used to simulate the model galaxies. 
The dark matter halo mass varies from 2 to 40 x Md 
which gives, on average, a mass resolution of ~ 10 5 M Q . 
We have observed that with such a mass resolution the 
shot noise in the simulation reduces considerably and 
bending instabilities do not grow sufficiently to induce 
an outer disk warp during the evolution. Each model 
galaxy in our simulation has be en evolved with the Gad- 
get code (Sp ringel et al. ||2001| ) which uses a variant of 
the leapfrog method for the time integration. The forces 
between the particles are calculated basically using the 
BH tree (with some modification) algorithm with a tol- 
erance parameter 9toi =0.7. The integration time step 
is ~ 0.4 Myr and each model galaxy is evolved for more 
than 5 Gyr in our simulation. The orbital time scale at 
the disk half mass radius is ~ 294 Myr for the model 
mk97 (see Table. [T]). For convenience, we describe the 
parameter space for the simulations below. 

3.1. Selection of the parameter space 

A wide range of parameters in the 3-D space spanned 
by [Q,o~ z /o~r,Mh/Md] is considered. The exact relation- 
ship between these parameters is dependent on how an 
active dark matter halo interacts with the other active 
disk and bulge components. However, we notice that 
if the central radial velocity dispersion is held constant, 
there appears to be a positive correlation between the 
Toomre Q parameter and the dark halo-to-disk mass ra- 
tio (Mh/Md) as expected from a simple understanding 
of the disk dynamics with a rigid dark halo. The initial 
thickness of the disks (h z /Rd) in our sample of model 
galaxies ranges from thin (h z /Rd ~ 0.1 or 0.15) to su- 
perthin {h z /Rd < 0.1). To give an example, the value 
of h z /R d ~ 0.07 in the superthin galaxy UGC 7321 (?). 
We explore a wide region of the parameter space spanned 
by Q, Mh/Md, cr z /a r (see Fig. [TJ and provide, here, the 
physical basis for our choice. Each model galaxy has a 
distinct Q and o~ z /a r radial profile; we quote the values 
of Q and a z /a r at the disk half mass radii (see Tabic [T]for 
representative models). Specifically, most of the galac- 
tic disks are vertically cold, but radially hot because we 
primarily aim to simulate low surface brightness (hence 
LSB) galaxies where the rotation curve is normally dom - 
inated by the dark matter halo (|de Blok et al. I [2001). 
In fact, many of the rotation curves in our sample (e.g., 
Fig. [2]) are such that the dark matter dominates the disk 
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Fig. 1. — The parameter space for all the models. Most of the 
galaxies in the sample are radially hot. 
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Fig. 2. — Rotation curve for the model galaxy mk97. The solid 
line is the total rotation curve, dashed line is due to the dark 
matter halo, dash-dot-dash line for the disk and dotted line is for 
the bulge. 

rotation curve from the center of the galaxy. It is also 
known that the stellar disks of LSB galaxies are thinner 
than thei r high surface brightness (hence HSB) coun- 
terparts (iBizvaev fc KaisinT I2004D . LSB galaxies are 
known to be poor in star formation activities despite the 
presence of neutral hydrogen gas as in otherwise nor- 
mal galaxies. This may indicate that these galaxies are 
dynamically hot (with high Q as most of the galaxies 
are in our sample) to prevent them from disk instabili- 
ties. Even if the disks suffer from instabilities, the result- 
ing non-axisymmetric features could be young. In fact, 
the low metallicities, high gas-to-star mass ratios, and 
blue colors of most LSB galaxies indicate that these sys- 
tems are probably youn ger than their HSB counterparts 
(|Vorobvovet al. II2009T) . 

4. HEATING MODEL 

We use the following simple expression (eq. 4) to fit 
all our simulation data, treating this functional form to 
quantify the effects of heating. It will be applied to the 
central region and to the outer region of the galaxy. 

<J z (r,t) = £r (r) + oi(r) x (t/T orb ) a . (4) 
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TABLE 1 

Initial conditions for the representative model galaxies 
(arranged in ascending order of Q) and the heating 

EXPONENTS COMPUTED AT THE END OF EACH SIMULATION. 



Galaxy 
models 


Q 




M h 




a ut 


mk51A 


1.21 


0.611 


6.52 


1.01 


0.326 


mk56A 


1.45 


0.594 


10.15 


2.59 


0.557 


mkl07 


1.66 


0.666 


9.83 


2.29 


2.00 


mkl04 


1.77 


1.480 


15.80 


1.67 


0.67 


mk97 


1.84 


0.447 


7.88 


2.31 


0.420 


mkl5 


2.06 


0.237 


5.34 


1.13 


0.809 


mk53 


2.31 


0.294 


11.97 


0.828 


0.289 


mk52 


2.33 


0.618 


14.47 


1.81 


0.707 


mkl7 


2.41 


0.151 


5.10 


0.724 


0.236 


dmA 


2.37 


0.433 


12.97 


0.776 


0.588 


mkll2 


2.47 


0.344 


9.44 


2.105 


0.215 


mk5 


2.50 


0.149 


9.85 


0.406 


0.448 


mk25 


2.53 


0.181 


7.07 


1.02 


0.459 


mk63 


2.64 


1.160 


14.12 


0.865 


0.902 


mk27 


2.71 


0.187 


7.16 


1.13 


0.30 


mk64 


3.00 


1.110 


16.50 


0.723 


1.39 


mk33 


3.01 


0.260 


8.75 


1.82 


0.41 


mk43 


3.05 


0.606 


11.52 


1.30 


1.63 


mk34 


3.15 


0.238 


6.50 


1.20 


0.597 


mk58 


3.16 


0.850 


16.75 


1.27 


1.02 


mk48 


3.33 


0.143 


10.47 


0.987 


0.27 


mk55 


3.73 


0.540 


27.60 


1.16 


1.10 



In the above equation, ctq and a\ are independent of time, 



and 



7~orb 



is the orbital time scale measured at the disk half 



mass radius, a is the heating exponent and is the loga- 
rithmic time derivative of the vertical velocity dispersion 
a = d\xia z /d\xit. We determine the time evolution of 
a z in the inner region (< lRd) and in the outer regime 
(~ 4i?d) of the disk. The definitions for the two regions of 
the disk are maintained throughout the paper. The heat- 
ing in each of these regions of the disk, characterized by 
apparently distinctive dynamical behaviors, is described 
by the resulting heating exponents given as and a ou t 
respectively. Higher values of a imply that the rate of 
heating is steeper as the galaxy evolves, whereas lower 
values indicate a gradual heating process. A thorough 
knowledge of these values are crucial for understanding 
the type of heating processes at work in the galaxy and 
possibly for identifying the mechanism responsible for the 
heating. The fitting is performed using a robust non- 
linear least squa res fitting algorith m due to Levenberg 
and Marquardt (jMarkwardt 1 120091: iMorel I1978D . The 
best fit model parameters are chosen based on the mini- 
mum value of the chi-square (Xmin)- Using the resulting 
parameters, an useful estimate for the time scale of heat- 
ing can be obtained from the above formula and this can 
also be compared with the relaxation time scale for the 
model under consideration. 

Let us consider that in time Th, the change in the verti- 
cal velocity dispersion Aa z ~ <7o (r) at a particular radius 
in the disk. Then differentiating eq.[4], it can be shown 
that 



Th 



ao\ (r) 



l/a 



7~orb • 



(5) 



We evaluate this time scale for different models pre- 
sented here. In comparison, the energy relaxation time 
scale for the N-body model with a total number of parti- 
cles Nj size R „ and softening e s can be written following 
iHunag et al~1 (|1993| ) as 



-Kl/2 

32^1og(i? s /e s ) 



Torbi 



(6) 



where we have used the crossing time scale T cross = 
/,! - T orbi assuming that the circular velocity V c re- 
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mains almost flat beyond R1/2', R1/2 denotes the half- 
mass radius. Considering N = 5 x 10 6 , R s = 10, 
R1/2 — 2.5 and softening parameter e s = 0.003, we have 

Trelax — 5.6 X 10^T or 5. 

Previously, disk heating has been modelled as a diffu- 
sion process in the velocity space. If the diffusion coeffi- 
cient D z remains constant in time, it c an be shown t hat 
the vertical velocity dispersion evolves (jLacev 1 11 9 9 ID as 



1/7 



f 

7~orb 



r 



(7) 



Note that the above formula is applicable to the radial 
dispersion also. The expon ent 7 = 1/2 denotes hea ting 
due to the halo black holes (|Lacey &: Ostrike"r"lll985l ) . In 
addition to eq. [4] , we also apply the diffusion model to 
fit our simulation data to some models for comparison 
and evaluate the corresponding diffusion coefficient D z . 
Let us consider that Td is the time scale over which the 
change in the vertical velocity dispersion Aa z ~ a z o- 
Then Td can be obtained as follows 



r d = (al /D z )T orb . 



(8) 



We denote Td as the diffusion time scale. In general, 
Td and r re i ax are comparable. Internal evolution of a 
disk galaxy is driven by various non-axisymmetric per- 
turbations such as bars or spirals which redistribute the 
energy and angular momentum within the disk and/or 
between the halo and the disk. The time scale for this 
process is longer than the dynamical time scale (Td yn ) 
and is usually known as the secular ev olution time scale 
(r sec ) (|Kormendv fc Kennicutt I [2004D . r sec could be a 
few to few 10s of orbital time scales (r or b). In ordering 
these time scales, T dyn < T orb < r sec < r d < T retax . We 
will use these definitions and compare the relevant time 
scales in the following to determine the relation between 
Th and the other time scales. 

5. NON-AXISYMMETRIC STRUCTURES AND DISK 
HEATING 

It is well established that stars undergo continuous 
heating in both the radial and vertical directions in galax- 
ies as they age. Although there has been substantial 
progress in understanding the radial heating, there is lit- 
tle progress for the vertical heating. Spiral structures 
have played a significant role in the studies of disk heat- 
ing in the radial direction. We note, however, that a 
steadily rotating quasi-st ationary spiral pattern heats the 
disk stars only negligibly (|Binnev fc Tremaine II 19871 ). In 
this case, sufficient heating of disk stars requires the 
spiral structures to be stochastic in n ature. Transi ent, 
swing-amplified spiral density waves (|Fuchs I l200l or 
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Fig. 3. — Particle image of mk97 showing the face-on and the 
two edge-on views of the initially axisymmetric disk (at T=0.) 

multiple spiral density waves (jMinchev fe Quillen l l2006) 
have recently been shown to be effective in the radial 
heating of disk stars. 

The primary mechanisms suggested for the vertical 
heating of the disk stars in early works range from a 
combin ation of stochastic spira ls and giant molecular 
clouds ([Jenkins fe Binnev1 ll990) to resonant heating of 
the trapped stars in the cha otic phase space in the pres- 
ence of a 3D bar potential (Pfennigcr 1985). More re- 
cently, other mechanisms have been suggested includ- 
ing a bending instability of the entire disk or of the 
bar (ISotnikova fe Rodionov 2003) and repea ted impact 
of globular clusters (jVande Putte et al. I [2009; , although 
it produces only a small change in the vertical velocity 
dispersion. Disk galaxies are rich in non-axisymmetric 
structures, which may form either by means of mergers, 
tidal interaction or internal instabilities. Examples in- 
clud e bars, transient sp irals, rings (e. g. NGC 1433, NGC 
6300 lButa etal~ll200l( ). lopsidedness (jSaha et al. Il200l . 
bending waves manifes ted in the form o f a warp (e.g. 
ES0121-G6, NGC 4565lSaha et al. ll2009h o r corrugation 
(e.g in IC 2233 iMatthews fe Uson 1 120081: iSaha et al. I 
l200a in NGC 5907) of the disk mid-plane. These non- 
axisymmetric instabilities often drive the disk away from 
equilibrium and are generally time dependent. Each of 
these features are, in principle, a potential source of 
heating of the disk stars. Since our understanding of 
the interactions of these features with the surrounding 
live dark matter halo is not well understood, it is diffi- 
cult to disentangle the signatures of each of these non- 
axisymmetric modes on the net heating of a group of 
stars. However, some progress has been made especially 
for our Galaxy in which it has been shown based on the 
analysis of the Hipparcos data that there exists a gap 
due to the 2 : 1 resonanc e with the bar in the Galactic 
UV plane (|Dehnen Il2000h . 

In the following, we study the nature of the vertical 



Fig. 4. — Particle image of mk97 at the end of the simulation 
(T=5.4 Gyr). The disk forms a peanut bulge which lasts till the 
end. 



0.75 - 




Fig. 5. — The amplitude of the bar is shown as a function of time. 
The solid line represents a type-I bar generated in the center of the 
model mk97. The dashed line indicates the amplitude of spirals in 
the outer region. 

heating of the disk stars due to bar growth and transient 
spirals. We also briefly discuss the effect of halo noise and 
the corrugation and/or bending waves on the net heating 
process. We carry out a series of N-body simulations to 
first systematically identify the sources of disk heating in 
our model galaxies which are free from tidal interaction. 
About 70 simulations have been performed covering a 
wide range of the above mentioned simulation space and 
a statistical attempt to characterize and understand the 
nature of vertical heating in various regions of the disk 
has been carried out. 

5.1. Vertical heating due to bar growth 

Although most of the galaxies in our sample are radi- 
ally hot, they form a bar at the center during its time 
evolution. Some of the model galaxies form a bar within 
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Fig. 6. — The vertical heating in the inner region (< li?^) of 
the galaxy model mk97 as a function of time. The solid line is the 
model fitted to the simulation data. cr z (0) is the initial value of 
the vertical dispersion at this radius. 
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Fig. 7. — Same as Figure [5] except for model mkf!2. 




Fig. 8. — Same as Figure [6] except for model mkfl2. 

about 3 to 4 rotation time scales, while for others it takes 
more than about 10 rotation time scales to reach nearly 
the same amplitude. It is found that a bar forms even 
for galaxies characterized by Q ~ 3.7 (model mk55) in 
our simulation, implying that it is difficult for a rotat- 
ing disk within a live dark matter halo to avoid forming 
bars at the center. Hence, bar formation appears to be 



quite common and generic in disk galaxies and a live 
dark matter h a lo pla ys a significant role. As shown by 
Athanassoula (20Q2|) a live dark matter halo supports, 
rather than suppresses, the bar instability in the disk. In 
contrast, bar formation in the disk could be suppressed in 
the presence of a non-responsive rigid dark matter halo 
(THohl II1976T) . The growth of the bar depends on the 
initial condition and grows via non-linear processes as 
the disk stars interact with the live dark matter halo 
particles. The nature of this interaction is yet to be 
complet ely understood a l though recent N-body simula- 
tions b v lDubinski et al. I fl2009ft . [Sellwood fc Debattista I 
(|2009D . and |Kly£in_eialJ ((20091) have revealed insights 
on the nature of bar dynamics in disk galaxies. For our 
choice of parameters and initial conditions, we do not 
find any unique trends for the bar growth in our sim- 
ulation. The sample of galaxies in our N-body simu- 
lations differs from typical galaxy models studied pre- 
viously (where Q is normally held constant throughout 
the disk or Q < 2.0 in general). Overall, we observe 
that the growth of bars in our sample falls under two 
broad distinct categories. In the first category, the bar 
grows quickly to a peak amplitude within ~ 5 rotation 
time scales which is about an order of magnitude less 
than r sec and nearly reaches saturation or starts grow- 
ing again (for convenience, we call these type- 1 bar). 
In the second category, the bar continues to grow but 
slowly as compared to type-I and does not show any ten- 
dency to saturate (we call these type-II bar). Type-II 
bars grow on a secular evolution time scale to reach the 
same amplitude as in type-I. These two main features of 
bar growth are illustrated in the appropriate figures be- 
low, although there are some intermediate cases present 
in our simulation. In many cases, the bar undergoes 
the w e ll known buckling ins t ability (iCombes fc Sanders 



19811: iPfenniger fc Friedhl 119911: iRaha et al. I 
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1991 



2006 



and references therein) following which the bar takes the 
form of a peanut and/or X shaped bulge. The buck- 
ling instability driven by the anisotropy in the veloc- 
ity dispersions is an important phase of the bar evo- 
lution as it leads to the form ation of a pseudo-bulge 
(|Kormendv fc Kennicutt I [2004D directly influencing the 
secular evolution of disk galaxies. However, the con- 
dition for the onset of this instability and the exact 
underlying cause of this phase remains a matter of 
further investigation. For example, recent work by 
IMartinez-Valpuesta fc Shlosman I ([20041 ) has shown that 
the buckling instability weakens the bar and the resulting 
peanut shaped phase lasts for several Gyrs in their simu- 
lation. We also find that the peanut shaped phase of the 
bar is long lasting, typically surviving till the end of our 
simulation (> 5 Gyr) producing boxy/peanut/X-shaped 
structure in the central region. The morphological evolu- 
tion of the disk in model mk97 is depicted through Fig. [3] 
and Fig. 2) These snapshots are taken at time T = and 
T = 5.4 Gyr which denote the beginning and the end of 
the simulation respectively. The initially axisymmetric 
disk evolves to form a bar which eventually transforms 
into a peanut morphology. Inspection of the evolution of 
the bar-amplitude reveals a diverse behavior while it is 
in the peanut phase; in some cases the amplitude con- 
tinually grows and in others, it saturates. While the 
bar strength continues to increase, the evolution of its 
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each N-body snapshot using the following formula: 



Fig. 9. — Solid line represents the amplitude of the bar (type-II) 
generated in the center of the model mkl07, while the dashed line 
is for the amplitude of the spiral in the outer region. 
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Fig. 10. — Vertical heating in the central region of the galaxy 
model mkl07. Solid line is the model fitted to the simulation data. 
0"z(O) is the initial value of the vertical dispersion at this radius. 



pattern speed is not always correlated in the sense that 
the pattern speed does not always decrease. In particu- 
lar, the bar pattern speed in some cases remains nearly 
unchanged with time while its amplitude grows (as in 
mkl7, mk33, mkl04). A similar anti-correlation between 
the bar gro wth and pattern speed evolu tion has been 
discussed by I Valenzuela fc Klypin (|2003l ) and more re- 
cently by iVilla- Vargas et al. I ( 2003) . However, there are 
also the 'normal' cases where the pattern speed decreases 
with time, while the bar grows to a higher amplitude, by 
losing angular momentum to the halo via dynamical fric- 
tion. We notice that the bar growth is slow (typically a 
factor of 5 to 10 less) in cases where the pattern speed 
is observed to remain nearly constant in time. Overall, 
the emerging diverse behaviour of the bar characteris- 
tics precludes detailed understanding for inferring gross 
physical properties (e.g. bar size, velocity dispersions) of 
the disk based on a few simulation results. As a conse- 
quence, the results of a large number of simulations have 
been used to approach the problem in a statistical sense. 

For each of the model galaxies, we compute the m = 2 
Fourier component (A 2 ) of the particle distribution from 



A, 



1 



(9) 



In the above equation 4>j(r, t) is the phase of the j th par- 
ticle at position r and time t and N r represents the num- 
ber of particles used. Since A m is a complex number, the 
modulus of the m th Fourier component is obtained by 

\A m \ = ^J$i(A m ) 2 + %(A m ) 2 (|Sellwood fc Debattista I 
2009) and the corresponding phase is given by (f> m = 
^ tan -1 3(A m )/5R(A m ), (m ^ 0). In general, the radial 
variation of A2 shows a pronounced peak corresponding 
to a bar in the central region of the disk and often a 
second peak indicating the presence of a spiral pertur- 
bation in the outer parts of the disk. The variation of 
the peak amplitude of A 2 in the inner region with time 
reveals the growth of a bar in the galaxy. The amplitude 
A 2 is normalized to the amplitude of the axisymmetric 
mode (m = Fourier component). In all the subsequent 
figures, we present the smoothed normalized A 2 . For a 
wide range of model bars, the following linear regression 
model is fit to the time evolution: 



log A 2 {r,t) = A° 2 +p\og{t/T orb ), 



(10) 



where the slope f3 is the logarithmic time derivative of 
the bar amplitude. The linear model in the log-log plot 
translates to the bar amplitude evolving as A 2 (r, t) — 
A%{t / T or bY . In the following, represents the growth 
rate of a bar and fi a represents the growth rate of a spiral 
(in the outer region). 

The time evolution of the normalized bar amplitude in 
model mk97 is illustrated in Fig. The bar amplitude 
nearly saturates through a series of high frequency oscil- 
lations while undergoing the buckling phase. Such oscil- 
lat ions have also been reported in previous simulations 
bv I Valenzuela fc Klypinl (|2003l ) and IVilla- Vargas et al. I 
(2009). Our estimate shows that the typical period of 
such oscillations is ~ 1/2 x r or b or less than T oro in 
general. This roughly implies that the amplitude oscil- 
lates with a frequency close to k, the radial epicyclic 
frequency. The vertical velocity dispersion (cr z ) of the 
stars in the inner region of the disk is plotted as a func- 
tion of time in Fig. |6] for model mk97. From Fig. [6l 
it can be seen that a z changes its slope at ~ 12 rota- 
tion time scales, nearly following the overall amplitude 
evolution of the bar. The best fit parameters of the 
heating model (eq.[4]) are a = 0.09, a± = 1.7 x 10~ 4 
and oiin = 2.31. Using these parameters in eq.(5), we 
obtain Th = 10.56 x T or \, <C T re i ax . The heating time 
scale due to the bar growth is found to lie in the range 
i~dyn < Th <C Treiax- We have also applied the diffusion 
model given by eq.(7) to the inner region of mk97 and 
find that the heating model given by eq.(4) provides a 
better fit than the diffusion model. In particular, the 
diffusion coefficient D z = 5.9 x 10~ 4 and 7 = 109 are 
found to be absurdly large. The resulting time scale 
Td = 1.5 X 10 3 T or6 . 

Interestingly, the temporal variation of a z does not re- 
flect any such oscillation while the radial velocity dis- 
persion undergoes oscillation with nearly the same fre- 
quency as the bar amplitude beyond about 10 rotation 
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time scales. 

The growth of a bar in model mkll2 is similar to that 
in model mk97, although the disk in model mkll2 is hot- 
ter (Q = 2.47) than mk97. Fig. [7] shows the normalized 
bar amplitude and the time evolution of a z in the disk 
inner region is shown in Fig. |8l The best-fit parameters 
from model mkll2 are found to be similar to mk97. 

In the case of model mkl07, the bar does not reach 
saturation and continues to grow as shown in Fig. [9l 
The vertical velocity dispersion in the inner region of 
the disk grows more gradually (Fig. [T0|) following the 
evolution of the bar. Note that the overall increase in 
the vertical velocity dispersion is less than that of model 
mk97 and mkll2. We find that the heating model given 
by eq.(4) better approximates the simulation data than 
the diffusion model. The resulting diffusion coefficient 
D z = 2.2 x 1CP 4 is similar to model mk97 and = 
4.4 x 10 3 T orb . For models mk97, mkl07 and mkll2, we 
find that, in general, the diffusion model overestimates 
the dispersion values in the inner region of the disk. 

In all our simulations, the vertical velocity dispersion 
increases following the growth of the bar. Note that all 
the three galaxy models (e.g. mk97, mkl07 and mkll2) 
undergo a phase of buckling instability during their evo- 
lution, but the disk stars are continually vertically heated 
prior to and during/after the buckling instability phase. 
This fact suggests that the vertical heating of the disk 
stars in our simulations is not strictly related to the bar 
buckling phase (e.g. the disk in model mkl7 does not go 
through buckling instability within the simulation time 
period but still shows continuous heating in the verti- 
cal direction). However, we note that the rate of vertical 
heating during the bar buckling phase is generally higher 
than in phases without such buckling. 

In the presence of a rotating bar (m = 2) potential, 
the planar motion of the disk stars can be coupled with 
the vertical oscillation (parametric resonance) at the lo- 
cation of the vertical resonances v/(Q,(r) — Qb) = rn/n, 
where v is the vertical frequency, fif, is the bar pattern 
speed, n is the number of vertical oscillations. For m = 2 
perturbation, the n = ±1 represent the vertical Lindblad 
resonances. For the n = ±2 resonances, the retrograde 
orbits in the inner region of the galaxy couple efficiently 
with th e vertical mot i on through t he 'B inney instability 
strips' (|Binnev Ill98lh . iPfenniger I (|1985t) has shown that 
indeed these 2 : 1 resonances inside the bar play a sig- 
nificant role in trapping stars, leading to rapid diffusion 
in the vertical direction. In the presence of a growing 
bar potential, the disk stars are subjected to a complex 
instability and the main effect of such an instability is 
to promote fast diffusion of the low angular momentum, 
high z amplitude and energetic stars to the nearby halo. 
The instability is shown to grow as the bar perturbation 
grows. Basically, the presence of a growing bar potential 
(i.e. a time dependent potential in the galaxy) breaks the 
time invariance symmetry (because bar growth is an irre- 
versible proprocess) and hence the jacobi integral which 
in turn enhances the ergodicity considerably a llowing the 
stars to vis it most of the chaotic phase space (|Pfenniger I 
119841 119851 ) . This provides a basic understanding of the 
heating mechanism in the presence of a growing bar. Of 
course, in N-body simulation the nature of bar growth is 
quite diverse as mentioned in the beginning of this sec- 
tion. The resonant heating of the disk stars in the pres- 
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Fig. 11. — The vertical heating in the outer region (~ &Rd) of 
the galaxy model mk97 as a function of time. The solid line is the 
same model fitted to the simulation data. a z (0) is the initial value 
of the vertical dispersion at this radius. 
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Fig. 12. — Same as figure ITT1 except for galaxy model mkll2. 

ence of a rotating (with changing pattern speed) growing 
(in amplitude) bar potential is probably more complex, 
leading to shifting locations of the resonances in the in- 
ner part of the galaxy and at the same time promising in 
the context of vertical heating in the disk. It is natural to 
examine how the vertical heating of the disk stars is cor- 
related with the growth of bar, the nature of the vertical 
heating due to a growing bar, and the typical exponent 
for the vertical heating. In section 6, we seek the pos- 
sible correlations between the bar growth rate (/if,) and 
the inner heating exponent (oj n ). 

5.2. Outer disk heating and transient spirals 

Strong two-armed spirals are difficult to excite and 
form in a radially hot disk galaxies. Using the Fourier 
decomposition of the disk surface density, however, many 
transient spirals (often with no steady pattern speed) are 
detected in our model galaxies. These transient spirals 
are often diffuse and generally weak in radially hot disks 
(Fig. E] and Fig. [7]) as compared to the relatively strong 
spiral formed in the radially cold disk model mkl07 (see 
Fig. [5]). These spirals are primarily confined to the outer 
region of the stellar disk, and we examine whether they 
play a role in the vertical heating of the disk stars in this 
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Fig. 13. — Vertical heating in the outer region (~ iRd) of the 
galaxy model mkl07. Solid line is the same model fitted to the 
simulation data. <r z (0) is the initial value of the vertical dispersion 
at this radius. 

region. In Figs . ITT | IT2l and IT3" 1 we present the time evolu- 
tion of the vertical velocity dispersion in the outer region 
of the disk for models mk97, mkll2 and mkl07 respec- 
tively. The slope of the heating curve in model mkl07 
(concave) is quite different from the other two models 
(convex) presented here. In the later stages of the evo- 
lution, <7 Z in model mkl07 increases faster than that in 
mk97 and mkll2 (where a z nearly saturates). Using the 
best-fit parameters of our heating model fitted to the 
outer region of mk97, we find that Th = 100t oi .& where 
cr = 0.025, ax = 8.6 x 10" 3 and a = 0.42. Applying the 
diffusion model in the outer region of mk97, we obtain a 
diffusion coefficient D z = 6.35 x 10 -8 and 7 = 0.214 and 
the corresponding time scale = 1.3 x 10 4 T or ;,, which 
is similar to the relaxation time scale (T re i ax ) quoted in 
section 4. In the case of model mkl07, the best-fit pa- 
rameters from our heating model yields = 33.7r or b- 
The diffusion model in this case gives D z = 8.48 x 10~ 4 
and Td = 1.15 x 10 3 r or b. These time scales indicate that 
the process of vertical heating in the outer region of the 
disk is slower in comparison to the inner region. We find 
that in the presence of a relatively strong spiral in the 
outer region, disk stars are heated to a greater degree 
in the vertical direction and the heating rate is found to 
be faster. The actual heating time scale, Th, (in Gyr) in 
mkl07 is lower by a factor of 6 in comparison to models 
mk97 or mkll2. Although, the exact mechanism through 
which stars are vertically heated in the outer region is not 
clear, growing spirals definitely play a significant role. 

The outer region of the stellar disk is not just simply 
described by the transient spirals, as corrugation waves 
or mild warps of the disk midplane are often present. 
It is known that the transient spirals can heat the disk 
efficiently in the radial direction, but poorly in the ver- 
tical direction. However, GMCs could redistribute the 
random energy of the stars in the vertical as well as in 
the radial direction through scattering processes. In the 
absence of GMCs or massive objects (such as halo black 
holes or dark clusters) in the outer parts of the galaxy, 
one seeks a process whereby the heating due to spirals or 
some other feature in the disk is redistributed in the ver- 
tical direction. Possible candidates, for example, include 
corrugation waves, large scale warps of the stellar disk 



or, in general, vertical motion of the disk stars coupled 
with the transient spirals. The growth and development 
of these bending waves or rather warps is known to a 
large extent to depend on the noise in the N-body sys- 
tem. We have verified that as the number of particles is 
increased (say by a factor of 10) in the system, the growth 
of these bending waves is significantly reduced because 
of the substantial decrement in the Poisson noise arising 
from particle discreteness. With an average mass reso- 
lution of < 10 5 Mq, only weak bending waves or mild 
warping of the stellar disk are present; comparatively 
strong warps can arise with a poorer mass resolution 
> 1O 6 M . Thus, it is unlikely that the bending waves 
alone can heat the disk stars in the outer region with 
higher mass resolution. However, the possibility exists 
that the non-linear coupling between the weak transient 
spira ls and the weak bending waves (jMasset fc Tagger"! 
1997) could redistribute energy in the vertical direction. 
As mentioned earlier, an interesting behaviour about the 
bars in radially hot disks inside a live dark matter halo 
was noticed. It was found that the type- 1 and type- II bars 
trigger transient spirals with distinctive characteristics in 
the outer disk. In the case of a type-II hai, the transient 
spiral continues to grow in amplitude, and this process 
continues as the bar grows (see Fig.[S]). This is found to 
be the case even in one of our hottest disk models (mk55 
for which Q = 3.73, see Table [T]). Because of the growing 
transient spirals, we find that the radial heating domi- 
nates over the vertical one throughout the disk. On the 
other hand, the amplitude of the spirals nearly saturates 
in the case of models mk97 and mkll2 (see Fig. [5] and 
Fig. [7|) where the bar is of type-I. Note that in this case, 
the radial heating nearly saturates throughout the disk 
(see Fig. [H] and Fig. [T5|) . Overall, there is a parallel de- 
velopment prevailing in the inner and the outer regions 
of the disk whereby a bar and spiral grow respectively. 

5.3. Anisotropic heating 

An important aspect associated with the heating of 
disk stars is its anisotropic nature as can be inferred 
from the fact that the observed ratio of the vertical to 
radi al velocity dispersion i.e. a z /a r < 1 (as in NGC 
488 lGerssen et al."1ll997l ). as found to occur for the so- 
lar neighbourhood in the Galaxy. This fact raises many 
questions. For example, why do disk galaxies have 
o~z/o~ r < 1? Can the heating mechanism preserve an 
initial a z /a r value? What causes the stars to be heated 
anisotropically? 

To gain some insight into these questions, we illustrate 
in Fig. Q3] and Fig.[TS]the evolution of the radial and ver- 
tical velocity dispersions in the inner and outer region of 
the model galaxy mkll2. It can be seen that the overall 
nature and the rate of heating is different for the radial 
and the vertical directions. In the inner region of mkll2, 
the vertical heating is more effective than the radial heat- 
ing and in the outer regime the radial heating saturates 
while the vertical heating continues slowly. We find this 
trend to be the case for most of the galaxy models which 
are radially hot and in which the bar evolves to satu- 
ration. On the other hand, for model mkl07 which is 
relatively cold (Q = 1.66), the radial heating dominates 
over the vertical heating throughout the disk due to the 
presence of a relatively strong spiral perturbation. It is 
clear that the non-axisymmetric structures (such as bars 
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Fig. 14. — Comparison of radial and vertical heating in the inner 
region in model mkll2. <r z (0) and <r r (0) are the values of initial 
dispersions at that radius. Solid line denotes vertical heating and 
dashed line the radial heating. 
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Fig. 15. — Comparison of radial and vertical heating in the outer 
region of model mkll2. tr z (0) and a>(0) are the values of initial 
dispersions at that radius. Solid line denotes vertical and dashed 
line the radial heating. 
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Fig. 16. — Vertical heating curves in the central region of the 
galaxies with different number of dark matter halo particle leading 
to different halo mass resolution. Solid line is for 0.5M, dashed line 
for 2M, dash-dot line for 3M, and dotted line for 5.0M particles 

or spirals) are essential for disk heating. 



5.4. Halo noise 

The role of a live or responsive d ark m atter halo has 
been emphasized bv I Athanassoula 1 ((2002) in the context 
of bar formation in N-body simulations. Unlike the case 
of a rigid dark matter halo, a live halo facilitates the 
growth of the bar in an N-body disk, likely through res- 
onant interactions with the disk particles. On the other 
hand, the finite number of particles employed in the N- 
body simulation can be an important source of disk heat- 
ing as well. Because of the fi nite number of p articles, 
two-body relaxation processes dVelazquez 1120051) can be 
effective in heating the disk stars efficiently. It is highly 
desirable to eliminate or, at least, minimize the Poisson 
fluctuation in the halo particle distribution by increasing 
the number of halo particles employed in the simulation 
as permitted by the available computational resources. 
To quantify this effect, we performed four additional sim- 
ulations on the galaxy model dmA (see Table [lj with dif- 
ferent number of halo particles (Nh). There is an order 
of magnitude variation in the number of halo particles 
between models with the poorest resolution (0.5 million 
(M) particles) and the highest resolution. The interme- 
diate two simulations use 2 and 3 million halo particles. 
In Fig. [TBI and Fig. [TT1 we show the four vertical heating 
curves in the inner region and in the outer region respec- 
tively. It is found that the vertical heating curves nearly 
converge for models with Nh > few x 10 6 . The conver- 
gence on the number of particles in our N-body simula- 
tion is in accor dance with the conclusi ons reached in the 
recent work by Dubins ki et al. I (|2009T ) in the context of 
bar formation and evolution. Comparison of Fig. 1161 and 
Fig. [17] reveals that the effect of halo noise is more im- 
portant in the central region than in the outer region. In 
order to quantify this trend, the time scales Th in the in- 
ner and the outer region were determined. We find that 
an order of magnitude increase in the halo particle num- 
ber increases Th by at least a factor of 10 in the inner 
region and ~ 2 in the outer region. In contrast, we find 
different behaviour when the temporal variation of the 
radial velocity dispersion in the disk was examined. Be- 
cause these models are dominated by spiral arms in the 
outer region and we find that by increasing the number 
of halo particles (from 0.5 M to 5 M) the spiral arms 
became even stronger. This resulted in a lower radial 
dispersion in models with lower number of halo particles 
compared to the models with higher number of halo par- 
ticles. In the case of model mkll2, the spiral arms are 
rather weak and the radial heating of the disk stars could 
have resulted from the halo noise. However, we find that 
the radial dispersion remains nearly constant over more 
than 10 rotation time scales (see Fig. H~4l and Fig. [15]). 
Based on these facts, it is possible to conclude that the 
noise due to the halo (with Nh > few million particles) 
is neither effective nor the dominant source of vertical 
heating in the outer region of the disk. 

In summary, we find a positive correlation between the 
growth of the bar at the galactic center and the vertical 
heating taking place in the central few kpc region. The 
formation of a bar is almost unavoidable in a stellar disk 
embedded in a live dark matter halo under the wide range 
of physical parameter space that has been explored. Bar 
formation occurs even in a very hot galaxy with its time 
for onset delayed at higher numerical resolution. In ad- 
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Fig. 17. — Vertical heating curves in the outer region (~ 4i?<j) 
of the galaxy models with different number of dark matter halo 
particles. Solid line is for 0.5M, dashed line for 2M, dash-dot line 
for 3M, and dotted line for 5.0M particles 
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Fig. 18. — Vertical heating exponent (ai n ) in the inner region 
(~ lRd) for all the models. Dependence on the bar growth rate 
(ft)- 



dition, analysis based on the four simulations on model 
dmA reveals that the strength of the spiral arms also in- 
creases with higher mass resolution. It has been noted 
that the vertical heating exponent is generally larger in 
the outer regions of a galaxy in the presence of strong 
spirals. Hence, the transient spirals in the 3D disk which 
are common in our simulations not only play a signifi- 
cant role in radial heating, but also the vertical heating 
of disk stars. 

6. CORRELATIONS 

In the absence of a clear understanding of the physi- 
cal processes responsible for the vertical heating of disk 
stars, we approach the problem by seeking possible cor- 
relations between the average physical properties (e.g., 
bar strength versus vertical heating rate) from our fairly 
large sample of N-body galaxies. To obtain a measure 
of the strength and slope of possible relations, we carry 
out a statistical correlation, measuring it between two 
variables x and y by the Pearson's product-moment cor- 
relation coefficient defined as follows: 



Corr(x, y) = 



cov(£, y) 
S x S y 



x)(yt - y) 



(N s - l)SJ y 



(11) 



In the above equation, cov represents the covariance 
and S denotes the standard deviation of the particular 
variable in the subscript. The summation is carried over 
the sample size N s . According to the Cauchy-Schwarz 
inequality the correlation can not exceed 1. If the value 
of the correlation is positive, it indicates an increasing 
linear relationship and in the case of negative value, it 
denotes a decreasing linear relationship. 

In Fig. [6j Fig. [U and Fig. 1101 the time evolution of the 
vertical velocity dispersion in the inner region of the disk 
is illustrated. In comparison with the evolution of the 
bar, it appears that the vertical velocity dispersion of the 
stars roughly follows the time evolution of the bar. Thus, 
for individual galaxies, there is a trend for the vertical 
heating, and it is interlinked with the growth of the bar. 
However, the exact relationship between the two and the 
reason for the vertical heating of the disk stars closely fol- 
lowing the growth of the bar is unknown. In the absence 
of a definitive answer, we seek to determine whether the 
trend found in individual galaxies is generic. For exam- 
ple, does it depend on the initial condition, and is this 
trend robust when examined on a large sample of galaxies 
with different initial conditions? We note that these are 
difficult questions to answer because galaxies with vary- 
ing initial con ditions evolve quite differently. As pointed 
out recently bv lSellwood fc Debattista I f|2009T ) a live dark 
matter halo interacts with the disk in a more complicated 
way than previously envisioned. 

In Fig-HH the dependence of the vertical heating expo- 
nent (ai n ) in the inner region of the galaxy on the growth 
rate of the bar ([3b) is illustrated for our sample of model 
N-body galaxies. The figure clearly shows that the heat- 
ing exponent in the inner region of the disk is strongly 
positively correlated with growth rate of the bar with a 
value of Corr(ai n , f3t) = 0.647. As pointed out earlier, 
when examined on a case by case basis, greater vertical 
heating results in galaxies which host stronger bars (e.g., 
models mk56A, mk52, mkll2). Our analysis suggests 
that as the strength of the bar increases, the disk stars 
continue to be heated vertically, irrespective of the bar's 
pattern speed evolution. The heating exponents tend to 
be higher in models with comparatively less massive dark 
matter halos. However, for more massive halos, the ver- 
tical heating exponent (oLi n ) tends close to unity. On the 
other hand, in the absence of a strong bar, the heating 
exponent is lower, cti n ~ 0.5. 

A weak correlation is found between the growth of spi- 
rals and the heating exponent in the outer region as 
can be seen from Fig. [19] corresponding to a value of 
Corr(a out , f3 s ) — 0.177. Its low value indicates that there 
may be more than one mechanism operating. Although 
the correlation between the heating exponents and the 
growth exponents for the spirals is not strong, higher 
exponents are associated with the presence of stronger 
spirals (e.g., model mk43, mk55, mk64) and the overall 
vertical heating is low when the spiral arms are more dif- 
fuse and weak (e.g., models mk5, mk51A, mk48). The 
vertical heating in the outer region is less extensive than 
the central region. The heating exponents in the outer 
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Fig. 19. — Vertical heating exponent (a «t) in the outer region 
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Fig. 20. — Correlation between the growth rates of bars (fit,) and 
spirals (/3 S ) for all the model galaxies. The figure shows that atleast 
weak spirals are formed in the disk whenever there is a bar. 

region in many models are very different from a ou t = 0.5 
indicating that the disk stars are not heated purely by 
massive objects in the halo or by the transient spirals, 
however, the exact source(s) which are responsible for 
the vertical heating in the outer region have yet to be 
identified. 

On the other hand, we find that there is a relatively 
strong correlation (see Fig. [20)) between the strength of 
the bars and the strength of the spirals in the disk with a 
value of Corr(f3b, f3 s ) = 0.565. This correlation indicates 
that in the presence of a bar there is always at least a 
weak spiral present in the outer disk. This may also 
suggest that such weak transient spirals are triggered by 
the presence of a bar in the disk central region. However, 
the vertical heating in the inner region and the outer 
region is not well correlated. As the correlation between 
the two is very weak Corr(a.i n ,a ou t) = 0.10, the source 
of vertical heating in the inner and the outer region of 
the disk is likely to differ. 

7. DISCUSSION 

We have studied the nature of vertical heating of the 
disk stars in a large sample of radially hot galaxies as a 
function of halo mass and o~ z / 'ay spanning a wide range. 
For an individual galaxy, we find that the stars in the in- 



ner region are heated vertically as long as a bar continues 
to grow. Such correlation is quite generic and found to be 
robust when examined on a larger sample of model galax- 
ies constructed from different initial conditions. Bars 
form in galaxies with diverse initial conditions but do 
not evolve in a comparable fashion likely reflecting their 
non-linear interaction with the surrounding dark matter 
distribution. However, we find that on a broad category 
there are two kinds of bar: type-I (growing rapidly and 
approaching saturation) and type-II (growing slowly with 
no evidence of saturation) as mentioned in section 5.1. 
Heating due to the growth of a bar differs from that due 
to a diffusion process since it occurs on a time scale of 
a few rotations in the disk central region. The diffusion 
time scale in a stellar system is similar to the relaxation 
time scale because diffusion mainly causes the stars ve- 
locity to change via two body encounters. Thus from our 
analysis, it can be said that the diffusion approximation 
is likely to be inappropriate for the vertical heating of 
disk stars in the central region. 

7.1. Growth of boxy/peanut bulge 

The boxy /peanut bulges are commonly believed 
to origin ate from the galactic bar 
heati ng dCombes et al. I Il990t 



19901: IPfennieer I 11991 



via the vertical 
IPfennieer fc Norman I 
lAthanassoulal 120051: 



Athanassoula fe Martinez- Valpuestal 120091 and ref- 
erences therein) and references therein). The growth 
of a bar in the N-body simulation with a live halo is 
a complex process because of the non-linearity in the 
growing process and its interaction with the dark matter 
halo particles. The growing bars in our simulation are, 
at times, associated with a nearly constant pattern speed 
and, at other times, with a decreasing pattern speed. 
Such a growing bar is likely to drive the inner disk to 
a non-equilibrium state and depending on its strength 
and pattern speed, the inner stellar disk may evolve 
on ~ a few dynamical time scales. Fundamentally, the 
growth of the bar facilitates the capturing of the disk 
stars into its various vertical resonances, e.g. the 1:1; 
2:1; 4:1 vertical ILRs (inner Lindblad resonances) etc. 
which are often pres ent in the bar region of the dis k 
(jCombes et al. |[l990t) . As shown by IPfennieer I (fl985h . 
the 4:1 vertical resonance becomes stronger as the bar 
grows, capturing disk stars to form the boxy bulge seen 
in many N body-simulations, although such studies do 
not explicitly include the interaction of the disk stars 
with a live dark matter halo. Furthermore, the presence 
of weakly dissipative orbits in a barred galaxy potential 
is shown to produce appreciable z-amplification as 
resonances a re crossed, hence facilitating the growth of 
boxy bulges (jPfenniger fc Norman 1 11 990). For example, 
a simple Ha miltonian mode l of the growth of a bar 
perturbation ([Quillen I 120021) illustrated the effect of 
resonant capture of the stars into the 2:1 vertical 
ILR, lifting them out of the disk plane to explain the 
growth of the boxy/peanut bulge. One common result 
amongst these studies is that the vertical resonances 
(when present) are efficient in capturing disk stars, 
facilitating their motion perpendicular to the disk plane, 
essential for the growth of a boxy bulge. We note that, 
in this context, the disk stars in our simulation are 
not only vertically heated at the resonance locations 
but throughout the bar region of the disk and often 
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with similar magnitude. This could imply either the 
resonant heating is not the only mechanism responsible 
for the overall vertical heating or broad resonances 
reflectin g strong chaos a re present throughout the bar 
region (jpfenniger I I1985T ). Although a growing bar 
with/ without a decreasing pattern speed could promote 
resonance sweeping through the bar region, it appears 
unlikely for the disk to accommodate the later scenario 
making it unclear how the non-resonant stars are 
vertically heated. In any case, it would be useful to have 
further insight into how the non-resonant stars might be 
contributing to the growth of the bulge. 

7.2. Superthin galaxies 

The disks of many galaxies in our sample are superthin 
and are dark matter dominated. From the galaxy forma- 
tion point of view, it is extremely important to under- 
stand how these superthin galaxies are formed and how 
they evolve. Can they preserve their initial superthin- 
ness? In the hierarchical structure formation scenario, 
a galaxy grows via a large number of mergers and as a 
result, the N-body simulations of galaxy formation un- 
derestimate the number of thin galaxies. So it remains 
a puzzle how superthin galaxies remain superthin. As- 
suming hydrostatic equilibrium, it can be shown that the 
disk half-thickness h z oc a z j \G p m id)l where p m id is the 
mid-plane volume density. Because of various unavoid- 
able heating processes o~ z would always increase with 
time and hence the disk thickness (h z ). Thus, either su- 
perthin galaxies somehow maintained their initial thick- 
ness or they evolved from an even thinner state to the 
present day superthin state. In the absence of environ- 
mental influences, the evolution of these galaxies would 
depend on the non-axisymmetric structures produced in 
the disk through internal instabilities. Our simulation 
shows that the non-axisymmetric structures such as bars 
or spirals are spontaneously formed in the disk and heat 
the disk stars in their respective regime of dominance. 
We briefly mention here that normally in radially hot 
and vertically very cold galaxies (e.g. mk25, mk27, mk48 
for which initial h z /Rd ~ 0.01 — 0.02), the disk thickness 
remains in the superthin regime with the final half thick- 
ness h z /Rd ~ 0.06 after 5 to 6 Gyr of evolution. Most 
of these galaxies form a type-II bar but there are ex- 
ceptions. If the initial thickness is greater than 0.03 or 
0.04, then the final thickness evolves to be in the 'thin' 
{h z /R d > 0.1) regime. 

7.3. Common heating mechanism? 

The vertical heating in the outer region is compara- 
tively more complicated because of the presence of spi- 
ral arms and, at times, mild warps or corrugations. The 
presence of these perturbations complicate the identifica- 
tion of the heating mechanism since more than one may 
operate simultaneously. It is known that a stationary spi- 
ral pattern does not lead to heating (except at resonance 
locations) and transient spirals heat more effectively in 
the radial direction. However, our analysis indicates that 
the transient spiral arms also grow (e.g., model mkl07) 



i.e. they are time dependent and, hence, a similar sit- 
uation may prevail in the outer disk as in the central 
region. In both regions, the disk stars respond to the 
time-dependent potentials. Thus, a common mechanism 
where the growing bars and spirals interact resonantly 
with the dark matter halo particles and contribute to the 
vertical heating in their respective regions of dominance 
may be operating. 

8. CONCLUSIONS 

Our investigation of the 70 N-body models of disk 
galaxies reveals a positive correlation between the growth 
of a bar and the vertical heating of the disk stars in the 
central region. A growing bar seems to contribute sig- 
nificantly to the vertical heating of the disk stars. Over- 
all, the heating exponent cti n > 1 for the various galaxy 
models studied. The disk stars in the central region are 
generally heated in the vertical direction by a factor ~ 3 
to 4 above their initial values over 5 to 6 Gyr. 

We find that the transient spirals are always present 
whenever a growing bar is present at the center. Most of 
these spirals are weak and diffuse, and it is likely due to 
the fact that the disks in our simulations are relatively 
radially hot. The numerical results show that the amount 
of vertical heating in the outer region is lower compared 
to the inner region with the disk stars in the outer region 
generally heated vertically by a factor of ~ 2 above their 
initial values over a time period of ~ 4 Gyr (if there is a 
growing spiral present) or more (otherwise). 

From the analysis of our simulations, we find that, in 
general, in radially hot galaxies the vertical heating in the 
central region dominates over radial heating and in the 
outer region the relative importance is reversed. In con- 
trast, radial heating dominates over the vertical heating 
throughout the disk for radially cold galaxies. We con- 
clude that heating due to non-axisymmetric structures 
appears to be most promising in the context of disk heat- 
ing problem in general. 

Our simulation results suggest that there is likely a 
common physical process through which the disk stars 
are heated vertically, which is active throughout the disk 
from inner to the outer region. Such a process should 
be investigated in detail by studying the vertical motion 
of the disk stars in the presence of a time dependent 
perturbing potential which could arise due to a growing 
bar or spiral arms inside a live dark matter halo. 

Acknowledgement 

All the simulation and analysis presented in this paper 
are carried out in 3 clusters namely the ASIAA computer 
cluster, TIARA cluster and the Academia Sinica Grid 
Computing (ASGC) center. The authors would like to 
express their thanks for the generous computing support 
from all of them. The authors also thank the referee, 
Daniel Pfenniger, for encouraging and useful comments 
which improved the presentation of the paper. KS would 
like to thank Yao-Huan Tseng, Sam Tseng, and Jason 
Shih for helping with the computer resources at various 
stages. KS would like to thank the Alexander von Hum- 
boldt foundation for its financial support during which 
the the paper was written. 



REFERENCES 



Ardi, E., Tsuchiya, T., & Burkcrt, A., 2003, ApJ, 596, 204 



Athanassoula, E., 2002, ApJ, 569, L83 



14 



Saha et al. 



Athanassoula, E., 2005, MNRAS, 358, 1477 
Athanassoula, E., & Martinez- Valpuesta, I. 2009, Chaos in 

Astronomy, Astrophysics and Space Science proceedings 

(Spinger-Verlag: Berlin), p. 77 
Barbanis, B. & Woltjer, L. 1967, ApJ, 150, 461 
Binney, J., 1981, MNRAS, 196, 455 

Binney, J., Dehnen, W., & Bertelli, G. 2000, MNRAS, 318, 658 
Binney, J. & Tremaine, S. 1987, Galactic Dynamics, Princeton 

Univ. Press. 
Bizyaev, D. & Kajsin, S. 2004, ApJ, 613, 886 
Buta, R. et al., 2001, AJ, 121, 225 
Carlberg, R. G., & Sellwood, J. A. 1985, ApJ, 292, 79 
Carr, B. J. & Lacey, C. G. 1987, ApJ, 316, 23 
Combes, F. & Sanders, R. H. 1981, A & A, 96, 164 
Combes, F., Debbasch, F., Friedli, D., & Pfenniger, D. 1990, A & 

A, 233, 82 
Debattista et al. 2006, ApJ, 645, 209 

de Blok, E., McGaugh, S., & Rubin, V. 2001, |astro-ph/01 07366 
Dehnen, W. 1998, AJ, 115, 238 
Dehnen, W. 2000, AJ, 119, 800 

Dubinski, J., Berentzen, I., & Shlosman, I. 2009, ApJ, 697, 293 
Edvardsson, B. et al. 1993, A & A, 275, 101 
Evans, N. W. 1993, MNRAS, 260, 191 

Font, A. S., Navarro, J. F., Stadel, J., & Quinn, T. 2001, ApJ, 
563, LI 

Fuchs, B. 2001, A & A, 368, 107 

Gerssen, J., Kuijken, K., & Merrifield, M. R., 1997, MNRAS, 288, 
618 

Gerssen, J., Kuijken, K., & Merrifield, M. R. 2000, MNRAS, 317, 
545 

Hanninen, J. & Flyn, C. 2002, MNRAS, 337, 731 
Hanninen, J. & Flyn, C. 2004, A & A, 421, 1001 
Hohl, F. 1976, AJ, 81, 30 

Hunag, S., Dubinski, J., & Carlberg, R. G. 1993, ApJ, 404, 73 
Jenkins, A. & Binney, J. 1990, MNRAS, 245, 305 
Klypin, A., Valenzuela, O., Colin, P. & Quinn, T. 2009, MNRAS, 
398, 1027 

Kormendy, J., & Kennicutt, R. C. Jr. 2004, ARA & A, 42, 603 
Kuijken, K. k, Dubinski, J. 1995, MNRAS, 277, 1341 
Lacey, C. G. 1984, MNRAS, 208, 687 
Lacey, C. G. & Ostriker, J. P. 1985, ApJ, 299, 633 
Lacey, C. G. 1991, in Dynamics of Disk Galaxies, ed. B. 
Sundelius, p. 257 



Matthews, L. D., & Uson, J. M., 2008, ApJ, 688, 237 
Markwardt, C. B. 2009, Proc. ADASS XVIII, ASP Conf. Ser., 

eds. D. Bohlender, P. Dowler & D. Durand, 411, 251 
Martinez- Valpuesta, I.& Shlosman, I., 2004, ApJ, 613, L29 
Martinez- Valpuesta, Shlosman, I., & Heller, C. 2006, ApJ, 637, 

214 

Masset, F., & Tagger, M. 1997, A& A, 318, 747 
Minchev, I. & Quillen, A. C. 2006, MNRAS, 368, 623 
More, J. 1978, "The Levenberg-Marquardt Algorithm: 

Implementation and Theory," in Numerical Analysis, ed. G. A. 

Watson (Springer- Ver lag: Berlin), 630, 105 
Pfenniger, D. 1984, A & A, 134, 373 
Pfenniger, D. 1985, A & A, 150, 112 
Pfenniger, D. & Friedli, D. 1991, A & A, 252, 75 
Pfenniger, D., & Norman, C. 1990, ApJ, 363, 391 
Pfenniger, D. 1993, in "Galactic bulges", IAU Symp. 153 (Herwig 

Dejonghe and Harm Jan Habing eds.) Kluwer, Dordrecht, p. 387 
Quillen, A. C. 2002, AJ, 124, 722 

Quinn, P. J., Hernquist, L., & Fullagar, D. P. 1993, ApJ, 403, 74 
Raha, N., Sellwood, J. A., James, R. A., & Kahn, F. D. 1991, 

Nature, 352,411 
Saha, K., Combes, F., & Jog, C.J., 2007, MNRAS, 382, 419 
Saha, K., de Jong, R., & Holwerda, B., 2009, MNRAS, 396, 409 
Sanchez-Salcedo, F. J. 1999, MNRAS, 303, 755 
Sellwood, J. A. & Debattista, V. P. 2009, MNRAS, 398, 1279 
Shapiro, K. L., Gerssen, J., & van der Marel, R. P. 2003, AJ, 126, 

2707 

Spitzer, L. Jr. & Schwarschild, M. 1951, ApJ, 114, 385 
Spitzer, L. Jr. & Schwarschild, M. 1953, ApJ, 118, 106 
Springel, V., Yoshida, N., & White, S. D. M. 2001, New A, 6, 79 
Sotnikova, N. Ya. &; Rodionov, S. A. 2003, AstL, 29, 321 
Vande Putte, D., Cropper, M., & Ferreras, I. 2009, MNRAS, 397, 
1587 

Valenzuela, O., & Klypin, A. 2003, MNRAS, 345, 406 
Velazquez, H. 2005, RMxAA, 41, 389 

Velazquez, H. & White, S. D. M. 1999, MNRAS, 304, 254 
Villa- Vargas, J., Shlosman, I., & Heller, C. 2009, ApJ, 707, 218 
Vorobyov, E. I., Shchekinov, Yu., Bizyaev, D., Bomans, D., & 
Dettmar, R.-J., 2009, A & A, 505, 483 



